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Abstract In this paper, a two-step regularization method is used to solve an 
ill-posed spherical pseudo-differential equation in the presence of noisy data. 
For the first step of regularization we approximate the data by means of a 
spherical polynomial that minimizes a functional with a penalty term consist¬ 
ing of the squared norm in a Sobolev space. The second step is a regularized 
collocation method. An error bound is obtained in the uniform norm, which 
is potentially smaller than that for either the noise reduction alone or the reg¬ 
ularized collocation alone. We discuss an a posteriori parameter choice, and 
present some numerical experiments, which support the claimed superiority of 
the two-step method. 
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1 Introduction 

Mathematical models appearing in geoscience commonly have the form of an 
ill-posed spherical pseudo-differential equation 

Ax = y, (1) 

where A is a pseudo-differential operator that relates continuous functions 
x £ C(Qr) and y £ C(Q P ) defined on concentric spheres Qr, Q p £ R 3 of radii 
R < p. For example, in satellite geodesy, this approach has been introduced 

(5j|25 


where the spheres C2r, Q p are models of respectively the surface of 


the Earth and the surface traversed by the satellite. 

Recall that a spherical pseudo-differential operator A : C{Qr) 
a linear operator that assigns to any x £ C(Or) a function 


C(f2 p ) is 


where 


oo 2fc+l ^ / \ 

Ax .— ^ ^ &k ^ ^ ( J ^ j 

k=0 j =1 ^ 


( 2 ) 


X k,j ~ \ R*k,j 




L2(^r) 


= rL, x 


( T ) Y k,j (^) 


are the spherical Fourier coefficients, and Ykj(-),j = 1,2,...,2 k+ 1, are the 
spherical harmonics 16 of degree k which are L 2 -orthonormal with respect to 
the unit sphere f2 1 6 R 3 , as a result of which 


( R Yk ’ j {r) ’ R Yk '’ j ' (r)) LAO r) ~ ( Yk 'i' Yk '^L 2 {n i) - SkvSjj'. (3) 

The sequence of real numbers («fc)fcLo referred to as the spherical symbol of 
A. We shall assume that a*, is positive, and converges monotonically to zero. 

In the case when the symbol sequence Ofe tends to zero fast enough the 
operator A is compact. Therefore, its inverse A -1 is unbounded, and keeping 
in mind Hadamard’s definition of a well-posed problem (existence, uniqueness, 
and continuity of inverse), we conclude that for this case the first and the third 
conditions are violated, and the problem 0-0 with y £ C{Rr) becomes ill- 
posed. Therefore, a regularization technique should be employed for solving 
it [|. 

As examples of the ill-posed problem 0, § we can mention the satellite- 
to-satellite tracking problem (SST-problem) with a*, = Ytl , the satellite 

gravity gradiometry problem (SGG-problem) with a*, = ( fc + 1 h fc+2 > ; etc. 

(for more details on these and other examples we can refer an interested reader 
to l>- These problems are severely ill-posed because of the occurrence of the 
geometric factor (R/p) k . 
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It is worth mentioning that many practical applications use a finite dimen¬ 
sional approximation of the solution o f (flj ), ©■ For example, Earth Gravity 
Models, such as EGM96 or EGM2008 [20| are parametrized by the spherical 
Fourier coefficients up to some prescribed degree M. 

Note that in applications the function y is assumed to be continuous. 
However, in practice one is provided just with a finite number of points 
{U}fL 1 C 17 P at which information about the values of y is collected. It should 
be noted also that the pointwise data y e (U) contain measurement errors, which 
can be modeled, for example, in the following way: 

I y e (ti) ~y(ti )| < e», i = l 

We assume for convenience that there exists a function y e £ C(fi p ) standing 
for the noisy version of the original function y, such that 


II V £ - y\\c(Q p ) < e := ^ax^ej, 

where are measurement errors. 

In such a setup the problem 0, § is reduced to the following spherical 
pseudo-differential operator equation 


M 


2/c+l 


AmX .— ^ ^ &k ^ ^ ^ k,j ( j — V • 

k =0 j =1 P 


( 4 ) 


Several regularization techniques can be used for treating For example, 
the method where the discretization level M plays the role of regularization 
parameter was discussed in |2j[3] 9 . However, in our case we assume that the 
value of M is prescribed, thus this approach cannot be used. 

The problem Q, but without noise in the right-hand side, was studied 


13 . In this noise-free case the solution x can be approximated by solving 


the equation 


A m x = V M y, 

where Vm ■ C(fi p ) —> P m{^ p ) is the so-called quasi-interpolatory operator. 
Here by Pjy(I2 p ) we denote the set of all spherical polynomials of degree less 
than or equal to M, or in other words the restriction to fi p of the polynomials 
in ]R 3 of degree less than or equal to M. Thus, in 


13 


the authors suggest 
constructing a polynomial approximation of y from the original pointwise data 
y(ti),i = 1,2 ,N, and then formally inverting the operator Am- 

In principle, this idea can be used also for the ill-posed case Q. However, 
in contrast to the approximation of a noise-free continuous function y on the 
sphere by means of polynomials, which has been discussed by many authors 
(see, for example, 7,24U26,28|30 ), the approximation of noisy functions y e has 


been studied only recently in |1[|21], to the best of our knowledge. Applying 
the method from [l 21 we first perform so-called data-smoothing (or noise 
reduction). After this data preprocessing the formal inversion of Am should 
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be safer. However, for performing the noise reduction step one needs a priori 
information about the smoothness of the function y, which is usually not 
available, a fact that makes the direct application of the scheme 21 not always 
appropriate. 

In a different direction, for estimating the Fourier coefficients Xkj directly 
from noisy measurements y e (U) a regularized collocation method has been re¬ 
cently presented in 18 . This method is based on the standard and widely used 
Tikhonov-Phillips regularization. However, it is well-known that the Tikhonov- 
Phillips method suffers from saturation 15 19 , meaning that the accuracy of 


reconstruction cannot be improved regardless of the smoothness of the solu¬ 
tion x. Another point is that while the Tikhonov-Phillips method has been 
well studied in the Hilbert space L 2 , to the best of our knowledge, no analysis 
has been done in the space of continuous functions, a natural choice for our 
noise model. 

In the present study we combine these two approaches. Moreover, in con¬ 
trast to the previous results, we will analyze the approximation in the space of 
continuous functions. A combination of two regularization methods into one 
can be seen as a two-step regularization, in which we use the composition 
Ra,M 0 T\,m of data smoothing operator Tx,m, and a regularized collocation 
operator R a ,M- In the literature there are not so many studies on two-step 
regularization, and we can refer only to [11||12| . However, the analysis in those 
papers does not correspond to the setting of our problem 0 @>. 

The paper is organized as follows. In the next section we present the regu¬ 
larization method for noise reduction and define the data smoothing operator 
T\ : m- In Section 3 we will give a short overview of the regularized collocation 
method, and define the operator R a ,M- Section 4 is devoted to theoretical error 
bounds for the constructed two-parameter regularization. We will show that 
the approximation has the potential to perform at least as well as the better 
of the one-parameter regularizations which are involved in the composition. 
Finally, in the last section we discuss an a posteriori parameter choice rule, 
and present some numerical experiments supporting the claimed superiority 
of our method. 


2 Data noise reduction 

At the first step of our scheme we approximate the noisy continuous function 
y e € C(f2 p ) by means of a spherical polynomial pm £ As discussed in 

the Introduction, instead of y we are provided only with pointwise measure¬ 
ments u e (tj),i = 1,2, Therefore we introduce the sampling operator 

S N :C{n p )\*2 for which 


W := {y e {ti),y e {t 2 ),---,y e {t N ))- 


By we denote the vector space equipped with the inner product 
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N 

i= 1 

1 /2 

the corresponding norm is ||??|| w = ( 77,77 )J . Here uji,u) 2 , ■ ■ ■ ,con are positive 
weights in a cubature formula, and t\, t^, ■ ■ ., tjv £ k? P are the corresponding 
cubature points, with the cubature rule having the property of being exact for 
all polynomials of degree up to 2M, 


N 

Vp £ P2 m(^ p ), E^) 

2=1 


[ p(Odn p (C). 

J Q p 


( 5 ) 


A method for approximating ?/ from S)v2/ £ £ based on such a cubature 
rule (but with all weights equal) was proposed recently in [l], while its gener¬ 


alization to other positive weights was presented in 21 
We briefly outline the method from 


21 and its extension to pseudo- 


differential operator equations. We start with the observation that the space 
Pjv/(f2p) of spherical polynomials p of degree at most M can be considered as 
the reproducing kernel Hilbert space (RKHS) Hk generated by the kernel 


M 


2k+l 


K{t,T) = Y J K 2 E 

k =0 j=l ^ 




t, T £ 12 p , 


( 6 ) 


where /3 = (/3 i,/?2, • • • ,/3m, • • ■) is a non-decreasing sequence of positive pa¬ 
rameters. Note that the inner product of H/c associated with this kernel can 
be written as 


M „ 2 2fc+l 

</^W = E§E(M-),/ 

fc=0 y 7=1 


i 2 (12 P ) 


Ykj - 


, 9 ) ■ (?) 

L*{n p ) 


The reproducing property (/, K(-, t))u k = /(r) for r £ fin can be verified 
easily. In this paper the numbers /?i, $2 > ■ ■ • > Pn are assumed to be given a 
priori. We shall see their role later. 

The approximation p m i n studied in [l 21 appears as the minimizer of the 
following functional 


Pmin = argmin|||5 Jv p-S'ivy e ||f, + A lbllE’ P £P m(^p)}> ( 8 ) 

where A > 0 is a regularization parameter. Note that © can be seen as the 
definition of the data smoothing operator T^m : C(12 p ) —> P m{^ p ) such that 
Pmin = T\, M y e - Note that the regularization term involves the 'H^-norm of p , 
and hence depends on the rate of growth of the parameters /3fc. The solution 
of ([8]) is given by the following theorem. 
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Theorem 1 Assume that the points {^} and weights {uji} are such that © 
holds. Then the minimizer T\,My e = Pmin in 0 has the form 


M 


2k+l 


w)(’) = E 


1 +A/3 1 

k—Q k '•= 


3 =1 


i—1 



(9) 


Proof It is known that the minimizer of the functional in ([8]) can be written 
as 


T x ,MV e = (A/ + S^Sn^S^Snv^ (10) 

where S% : -» TL K is the adjoint of S N , and I is the identity operator in 

TLk- By definition, for £ Rff we have 

N N 

(v, SNf)„ = = J2vi (K{U, ■), f(-))u K Ui 

i =1 2=1 

= ( '52 UJ i K ( t ir)Vi,f(-) \ 

\*= 1 / Hk 

thus there holds 


N 

2=1 

and thereby 


N 

{s* N s N f){-) = Y,“iK{t u -)m. 

2=1 

Inserting this expression into ( [To] ) , we observe that p m i n solves 


N N 

^Pmin(') T ^ ^Pmin^) = kJjK(tj , •)?/ (tj) (H) 

2=1 2=1 

The spherical harmonic expansion of the polynomial p m i n = T\,MV e is 


M 2/c+l / ^ 

Pmin(^) = ^ ^ ^ ^ \Pmim ~^k,j 
k—0 j =1 ' ^ 



) L 2 (r2 P ) 




t £ 17 P . (12) 


To find the coefficients in this expansion, we write the second term on the 
left-hand side of (111, on using ([5]) and then ([fi]), as 


N 


y^Wi/\(t;,t)p m in(ii) — (pmin, K(-,t)) 

L 2 {Q P ) 


2=1 


N 2fc+l 


= E^ 2 E (Pmin, Y k , 


k=o i=l 
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After expanding the right-hand side of (111 using (|g]), and then equating the 


coefficients of 


Yk 'i (p)’ 


^ , we obtain 


(A + /3 fc ) ( Pmin, Ifej ( p 




N ft\ 

Bk 'y ' WjYk,j ( — J y (t?;), t £ fl p 

»=l 'P/ 


which on substituting into (121 yields the desired result. 


At this point the solution of the first step depends on the regularization 
parameter A, and the penalization weights Bk- As mentioned in the Introduc¬ 
tion, the choice of the regularization parameter A will be addressed in the last 
section. A data-driven choice of the penalization weights /3k has been recently 
discussed in [2l] . 

The assumption that the function x on 12 r is continuous implies that 
x £ L 2 (1?r), and hence that its Fourier coefficients (jjYk,j (r) , x ) l 2 (o r ) w ith 
respect to the basis of spherical harmonics are square-summable, i.e, 


oo 2/c+l 

EE 

k —0 j =1 






2 

< OO. 


Any additional smoothness of x can be measured in terms of the summability 
of Fourier coefficients with some increasing weights depending on the sequence 
(Bk) or, as it is usual for the regularization theory (see, e. g., [14]), on the 
symbol (a*,). Therefore, it is convenient to introduce two Sobolev spaces W^, 
and W^ ,a , the first depending on the sequence (Bk) and the second on the 
symbol (a*-), and defined by 


6 3 ) 2 (« 5 ‘j(fl) lS )l2(BH) 

w +* := g e L\Qr) : 22 E — 

k =0 3 =1 




\\w6,3 < oo > , (13) 


I „ (k) ’ 9 ) l (n J 

W+’ a := { 9 6 L\Qr) : E ~ 

k =0 j = 1 


0 2 (®k) 


II W"^ 


,.<«>), (14) 


where 0,0 are non-decreasing functions such that 0(0) = 0 and 0(0) = 
0. In the literature, see, e.g., 14 , the functions 0,0 go under the name of 
“smoothness index functions”, and are usually unknown. 


3 Regularized collocation method after noise reduction 

After the first step of our method the reduced original equation Q is trans¬ 
formed into the equation 
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A m x = T x ,My e - (15) 

The regularized collocation method (see [18]) is now applied to this equation, 
yielding an approximate solution x e a x of the equation (151, defined as the 
minimizer of the functional 


\\Bn,mx - SNT\ t My e 

where Bn,m '■= SnAm ■ L 2 (Hr) 
parameter for the second step. The minimizer of (161 can be written in the 
form 


e, 


^ a \\ x h 2 (n R ), ( 16 ) 

and a > 0 is the regularization 


where B 


N,M ■ 


X a,\ ~ i a ^ + B n m Bj^ m) 1 B N M S]\[T\ My e , 

L 2 (Hr) is the adjoint of .B/v.m, given by 


M 2fe+l . N i 

(Bn,mV)(') = a k E] (^) 


fc =0 


3 = 1 


*= 1 


from which it follows, using Q, that 


M 2fc+l ^ 

n,mBn,mx = 'y ) a k y ] (^>) x k,j- 


B„, m ...... 

k—0 j =1 

Using (|l7| and ^ we then obtain explicitly 


m, 


(17) 


M 


u : = E 


dfe 


a + 1 + A/3p 


2fc+i 1 iv .. , 

E R Y ^ ( r ) 


fc=0 ~ ' ~*fc j=1 - —■ i=1 

If we now define R a ,M ■ Pm{^ p ) —> P m(^r) as 

(R a ,MP) (') := (( a/ + Bn 9 m b n,m) (•) 

M 2fc+l 1 AT 

a + a 2 (/?) E 


= E 

k—0 

M 

= E 


2/ e (^)- (18) 




AC j = l 

2fc+l 


i—1 


,j 

P 


k—0 


a + a 


k a — 


3=1 


- s y * J W ^ y * J ’ p , LAn ,, 


P(*i) (19) 

( 20 ) 


then we can write A = R a ,MT x ^My <L ■ Thus 7 ? Qj m reflects the regularized 
collocation second step of the method, whereas T XM corresponds to the noise 
reduction first step. 

As one can see from (181, in the end we approximate the unknown solution 
x on Hr by means of a spherical polynomial of degree M. In this context a 
question arises about the best approximation of a continuous function x by 
means of a polynomial of degree M. In turn the quality of best polynomial 
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approximation is determined by the smoothness of x. We shall assume that x 
lies in the intersection of and W^' a , see (131 or (14), thus the smoothness 


of x is encoded in <j> and /3 k, or i/> and a*,. For example, if the smoothness index 
function </>{t) and the sequence (3 = {/3k} increase polynomially with t and k, 
then Jackson’s theorem on the sphere (see j22], Theorem 3.3) tells us that for 
x £ there is v > 0 such that 


JgJ x - P\\ C (a R ) = O (M "). 


( 21 ) 


At the same time, if the sequence /3 = {/3k} increases exponentially then 
for polynomially increasing (/> and x £ we have 

™ ? i J x -P\\c(n R )=0(e-« M ), 

where q is some positive number that does not depend on M. 

In the error analysis of the next section we shall make use of a constructive 


polynomial approximation introduced in 26 , in which x £ C(f2u) is approxi¬ 
mated by Vmx £ P m(^r), given by 


M / , \ 2/c+l / f \ / 1 \ 

(Vmx) (f) = S'* (m) S e Y ^ (h) (r) ’ tel2x ’ 

( 22 ) 

where h (a “filter function”) is a continuously differentiable function on R + 
satisfying 

h(t) = { l’ \ 6 [?’ 1// ? ’ 0 < h(t) < 1 fort e R+ 

i 0, t £ [l,oo), 

Explicit examples of suitable filter functions h can be found in 271. It is im¬ 


portant for our later analysis that, as shown in 26 , Vmx is (up to a constant) 
an optimal approximation in the uniform norm, in the sense that, 


\\x-V M x\\ c{Qr) < c p jf M/2] W x -pWc(q r ) - 


(23) 


where c is a generic constant, which may take different values at different 
occurrences, and [•] denotes the floor function. 

In view of (21), for polynomially increasing </>, (3 and for x £ we have 


||x - V M x\\c{a R ) ^ c [ M / 2 ] " < cM "■ 

On the other hand, for exponentially increasing /3 and polynomially increasing 


(f> the theory 23 suggests taking h(t) = 1 ,t £ [0,1]. In this case the right-hand 
side of (231 has to be modified by multiplying by \[M, and by replacing [M/2] 
by M, thus for x £ there holds 


\x-V m x\\ c{Dr) < cSMJg" \\x-p\\ c{Qr) < cVMe qM . 
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4 Error estimation in uniform norm 


In this section we will estimate the uniform error of approximation of x by the 


polynomial x e a x given by (181. It is clear that 


c( n R) = \\ x - R ‘*,M T XM ye \\ c( n R ) 


(24) 


< \\ x ~ Vm x \\c(s 7 r ) 4 \\V M x ~ R a ,MU M y\\ C ( nR) 

+ ||(-Rq;,M — Ra,MTx t M)UMU 
+ t?CK,MTx,M 


'C{Q p )—^C{Qr) 


\C(Q r ) 

(ll Umv - y\\c(n p ) + II v - y e \\c(n p )) > 


where UmU is the same as VmX in (221 except that it is an approximation on 
the larger sphere 12 P instead of on 


M 

{U M y) (i) = £ 

k—Q 


7 \ 2fc+l 1 / 

- 

jri p \p 


Y kd (- ,»(■)> • (25) 

PJ / L2(l? p ) 


It is natural to assume that \\UmU — y\\c(n ) < e > since otherwise data noise 
is dominated by the approximation error and no regularization is required. 
We also restrict ourselves to the case when ||a; — Vmx\\ C (q r -) < ce, otherwise 
the term ||a; — VMX^c(n R )i representing the error of almost best polynomial 
approximation, will dominate the error bound. 


Then the bound (241 can be reduced to the following one 


||z - ®a,A|| C {O r ) ~ W VmX ~ R a,M U My\\ C (n R ) ( 26 ) 

+ II ( R a,M ~ R a,MT\ t M)UMy\\c(n R ) d~ C (1 "b W R a,MT\ t M\\c(n p )^>C(0 R )') e ' 

An estimate of the coefficient in the last term is given by the following 
theorem. 


Theorem 2 Under the conditions of Theorem 1 

N M 

yy (2 k + 1 )a k 


VMM C( 0 P )^C( 0 R ) i Rp ££_ 


' 4jr(a + a?)( 1 + A/3?) 


Pk 




i =1 fc=0 

where P k are the Legendre polynomials of degree k. 

Proof In view of the definition © and ( |T9) , together with ([3]) with R replaced 
by p and the property <§ of the cubature rule, we can write 

(R a , M Tx, M y) (t) 

M 


2fc+l 




k—0 


M 


1 v 2k + 1 elk 
k =0 


R p 4n a + a k 1 + A p k — 


N ft-t\ 

Yl UiPk e 
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where in the last step we used the addition theorem for the spherical harmonics 

ED 

2fc+l 


E Y kd( x ) Y k,j(y) = 2k ^ 1 Pk(x ■ y), x, y G Qi. 
i=i 7r 


Finally, we estimate the first two terms of the error bound (26). 


Theorem 3 Assume that the smoothness index functions <j>, ip are such that 
the functions t/<p(t), and t/ip(t) are monotone. Then for x € W^’P P| W^’’ a 
and Ax = y we have 


||VivfX Ra,MUfriy || A ||(-Ro,M Ra,AlT\^ Al) UMy ||) 

< cM[ip(a) \\x\\ w ^, a + 4>{\) ||x||^], 

where </>(A) = </>(A ),ip{a) = ip{a), ift/cppt) and t/ip(t) are non-decreasing, and 
4>(\) = X 1 ip(a) = a otherwise. 

Proof It follows easily from Q and ^ that 


ak \R 


YkJ (r)’ x 


L>2 (^r) 


= - Y *d k M 


L*2 (ftp) 


Hence from (22) and (25) we obtain 

M ' k \ ii 


(V M X) (t) = Y, h ( M ) ^ E ( R ) ( ~ n Y ^ ( n ) 


k—0 
M 1 2fc+l 


M) ak ' R 

J= i 


L 2 (r? p ) 




2 nj 


, ( UmV ) 


L 2 (O p ) 


Moreover, in view of (261, (251 and (|5| we also have the representations 

M 2fc+l N 

(RcmUmv) (*) = Y1 E lf Yk ’i ( p) (ti) 


k =0 
M 


J=1 

2fc+l 


I/ 2 (i7 p ) 


-E^E (s) (p) 'W» 

k=0 j=l 

while from ([7]) and ([3]) and two uses of Q we obtain 

M 2fe+l N 

(R„,MT x , M U M y)(t) = EE^rrW E ^ (£) E“*^ (5) 

k =0 j = l i=l 

M 2fe + l 

- E E (I) (h“ 


n 2 (n„) 
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Then we can write, using again (25), 


||VjvfX R a ,MU M y\\c^ E 


M 


E 


ak 


k =o 
M 

Y, h 

k =0 


a k a + ai 


2fc+l 


k / a— 


3 =1 


E jjM 


2fc+l 


a fc (a + aj;) ^ 


RJ \p 
t \ / 1, 


RJ \p 


-Jkj 


p)' v ) 


L^O p ) 


L 2 (n p ) 


c(n R ) 


c(n R ) 


Now using the property that 0</i(^)<l and then the Nikolskii in¬ 
equality (see, e.g., 17 , Proposition 2.5), we obtain, for x € W^' a , 


< 


II Vmx — R a ,M^My\\ c ^ nR ^ 


M 2fc+l 1 

e^eK 


k—0 


* .7 — 1 


i? 


R \R 


Y k 


4 «) •*) 


L 2 (O r ) 


c(n R ) 


< cM 


= cM 


M 2k-\-l 1 / , \ / i \ 

E^E^(i)(±Ms ,*) 

fc=0 k 3=1 \ / \ / 


fc =0 

M 




i 2 (« R ) 


Y- a ., 2 ^ 2 \~' 1 (k) K ) 

E^«)E- 

/ c —0 K .7 — 1 


4K) 


L 2 (O r ) 


< cM sup 


a 


a + u 


4(4) 


Fll w** 


The proof can be completed by using Proposition 2.7 from 14 . For con¬ 
venience, we present the corresponding argument here. 

Let us consider two cases. First, we assume that t/il>(t) is a non-decreasing 
function. Then keeping in mind that the function ^ is increasing, for 0 < u < a 
we obtain 


sup 

u<a Oi + u 

If a < u we can write 


4(4) 


< sup ip( u ) < 4( a )- 

u<a. 


sup 

cx.<.u 


a + u 


4(4 


< sup 

a<u 

< 


inf Q < 


au 4( u ) 
a + u u 
a 

- < 


< sup 

a<n 

a 


au 

a + u 


= ip(a). 


— u ip(u) ^(ai) 


sup 

OL<.U 


4(4) 

U 
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Assume now that t/^){t) is a decreasing function. Then 


sup 

u<af 


iiu) 


< sup 

u<af 


au ip(u) 
a + u u 


a + u 

Combining everything together we obtain 


^ V’K) 

< a — 7 \— < col. 
ai 


\\V M x - R a , M U M y\\ c{OR) < cM'tp{a) ||a;|| wV ,, 0 . 

The second part of the regularization error can be bounded by applying 
similar steps and using the assumption of the theorem that x £ . □ 


From (26 ) and Theorems 2 and 3 we can estimate the error of the approx¬ 


imation of x with the use of the proposed two-step regularization scheme as 
follows: 


I* - xC a ,\\\c(n R ) - ciW M A ) \\ x \\w^f> + $( a ) IMIwv.o] 

(2 k + 1 )afc 1 


(27) 


c max 


M 

E 


te ° R S, 47r (°fc + «)(! + A /3fc) Rp 


N 


2> A (w 

2=1 N 


e. 


From the error decomposition (27 ) one can easily obtain the error bounds 


for the single parameter regularization schemes involved in the two-step com¬ 
bination. If A = 0, which means that no noise reduction has been done, we 
derive the following error bound in the uniform norm: 


(28) 


On the other hand, if a = 0, which means that direct inversion after data 
pre-smoothing has been done, then 

(29) 


Now we illustrate the potential advantage of the two-step approximation 
compared to the single-parameter schemes. Assume that <j>(t) = t^, and ip(t) = 
t v for some fi, v £ (0,1). In this case the error bound ( [27] ) is reduced to the 
following one 


\ x ~ x a,x\\c(n R ) ^ cM< K A ) \\ x \\ww 

2k+ 1 


M 

c max Y 
teO R f—' 47rafc(l 


k —0 


Wl) Rp 


N 

^ ' LdiPk 
i=l 


t-tj \ 

Rp) 


:,A || c(o R ) - cM ^( a ) INI 

(2 k + l)a k 1 


c max 
t&n R 


M 

E 

fc=o 


(a + ay Rp 


N 

E 


'X i Rk 


f t-tj 

V Rp 
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X - x a,x\\ C {0 R ) ~ cM ( A ^ ll X ll W*-l> + IMIv^.») 


M 

+ c max > 


(2k + 1 )afe 1 

(a + a|)(l + A/3|) Rp 


N 

^ ' kJiPk 
i =1 



(30) 


Consider for simplicity the situation in which for fixed p and v, a = \ ,J ' V . 
Then the first terms of (28), (291 and (301 are essentially the same, whereas 
the last term in (301 is obviously smaller than that in either (28) or (291. 


5 Numerical illustrations 


In all our experiments the regularization parameters A, a are chosen by means 
of the so-called quasi-optimality criterion. This heuristic approach was origi¬ 


nally proposed in 29 and has been recently advocated in 10 


In the case of the noise reduction method <[8j), the quasi-optimality criterion 
chooses a regularization parameter A = A i from a set 


Qx, = {A = Ai = A 0 q\ * = 0,1,2,..., Jb} , q > 1, 


(31) 


such that 


II A m( t >-iMV' l ~ T Xl _ lMV *)\\ = min j 11(Tx i ,mv‘ ~ T\ i _ 1 ,M v‘ ) 11 , i = 0,1,. .. , i} , 

where ||-|| denotes the uniform norm in the space C = C(fln) of continuous 
functions. 

In the similar way we can apply the quasi-optimality criterion in the regu¬ 
larized collocation method. Then a regularization parameter a = is chosen 
from a set 


Qk = { a = a j = a or J , j = 0 , 1 , 2 r> 1 , 

in such a way that 


(32) 


II Ra k ,My e - -Rat- 1 , MU e || = min {11 R aj , M y e - R aj ^ lt My e \\ , j = 0, 1, ... ,K} . 

For the two-step regularization the quasi-optimality criterion can be im¬ 
plemented as follows. First, for every a = aj £ Q° K we choose A = A; = A(<x,) 
from the set (31) such that 




. r 

_ £ „e 

min < 

^QLj,\i ^aj ,Ai_i 



Next, we apply the quasi-optimality criterion to the sequence 
parametrized by aj £ Q More specifically, we select ak £ Q^ such that 
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xl i\ — x*L \ t „ \ — min -I x* \ — x* , j = 0, 1, . . . , K > . 

Then, a regularized approximate solution x e a A of our choice is defined by ( |T8j | 
with A = A(afe) and a = a*,. 

and assume that is the set of 

Gauss-Legendre points, for which the positive quadrature weights are known 
analytically. In this case N = 2(M + l) 2 , and we take M = 30. 

The data are simulated in the following way. First we generate a spherical 
function 


In all our experiments we follow 6,18 


M 2k+l 

y(t ) = ak 

k —0 j= 1 

where Xkj = (k + l/2)~ v gkj, k = 0 j = l,...2fc + l, v > 0, and 
gkj are uniformly distributed random values from [—1,1], This means that 
the solutions of the simulated problems have Sobolev smoothness v, and y 
can be seen as a result of the action of the operator A on a function x from 
the spherical Sobolev space H%. A noisy spherical function y e is simulated by 
adding a Gaussian white noise of intensity e = 0.05 to the values of the initial 
function y at the points 

To assess the performance of the considered schemes we measure the rela¬ 
tive error 


( p ) > * e 


where 


M 2k+l 


^r y ^(r)’ 


k =0 j—1 

and the approximation x e a x is given by |l8| . 

The results are displayed in Figure [Tj where each subfigure corresponds 
to different values of ak,v, and penalization weights /3k- In Figure [T] each 
circle exhibits a value of the relative error in solving the problem with one of 
10 simulated data, for each of three methods: regularized collocation method 
corresponding to the case when A = 0, noise reduction method with direct 
inversion (the case when a = 0), and the two-step regularization method. Note 
that such a form of graphical representation of the performance of different 
regularization algorithms is rather common (see, e.g., HP- 

In our experiments we tried to cover different degrees of ill-posedness of 0. 
([2]), i.e. the so-called moderate ill-posedness when the spherical symbol tends 
to zero polynomially (Figures [lj, [lji, [T ;), and severely ill-posed case with 
exponential symbol decrease (Figures |TJl, ^). Note that the spherical pseudo¬ 
differential operators with polynomial symbol have been studied in 13 , while 
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(a) Relative errors for the case a & 
(1.48) _fe , a = 3/2,/3 2 = a~\k 
1,2,..., M. 


(b) Relative errors for the case a^ = 
(1.48)-*> = 11/2, 0l = a~\k + 
1/2 ) 7/2 ,fc = 1,2, ...,M. 


ct=0i- OOOSIOO o 


A=0 • O o o 


O CD O O 00 


O 008800® O 


0.01 0.02 0.03 0.04 0.05 0.06 


1.5 2 2.5 3 3.5 4 4.5 5 


(c) Relative errors for the case a k = (d) Relative errors for the case a k = 
(k + l)~ 2 ,v = 3/2,/3 2 = a^ 1 , k = (k+l)~ 2 ,v = 11/2 , p 2 = a fe 1 (fc + 
1, 2,..., M. l/2) 7 / 2 ,fc = 1,2, 


o oo oo 
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(e) Relative errors for the case a & = 
(k + l) _2 ,t; = 11/2,/3| = a~ 1 (k + 
1/2) 11 / 2 , k = 1,2 


Fig. 1: Numerical illustrations. The comparison of the performance of single¬ 
parameter regularizations and their combination (Multi). In each of the figures 
the results corresponding to the two-step regularization are displayed at the 
top, the case of noise reduction with direct inversion is shown in the middle, 
and the case of regularized collocation corresponds to the bottom row. The 
horizontal axes show values of the relative errors (the vertical axes have no 
significance) 
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the operators with exponentialy decreasing symbol, appearing in geoscience, 
have been discussed in the Introduction. 

It is instructive to look at Figures [lji, [l]i and from the view point 

of 21 , where a priori choice of the penalization weights /3 k has been dis¬ 


cussed. For the situations described in the above mentioned figures the choice 
of /3fc suggested in 21 would be f3\ = a^ 1 (k + 1)". But such a choice re¬ 


quires knowledge of the smoothness v of the unknown solutions. In practice 
the smooothness v is usually unknown, and the Figures [l]i [lj:l correspond to 
the case of “underpenalization”. 

As can be seen from Figures MP in such a case it is not clear a priori 
which of the considered single-parameter regularizations performs better. At 
the same time, these figures show that the proposed two-step method automat¬ 
ically follows the leader. The situation does not change when the penalization 
weights are chosen according to [21] (see Figure [T^). 

In all our experiments the choice of the regularization parameters was 
done by using the quasi-optimality criterion described above with ao = Ao = 
1.78e — 5, q = r = 1.25, and L = K = 50. Note also that we add the term 
A = a = 0 to each of the sets (311, and (321. 

These numerical results confirm the theoretical conclusion that the con¬ 
structed two-step regularization methods can perform similarly to (or even 
better than) the best of single-parameter schemes. 


6 Conclusion 


We studied a two-parameter regularization scheme, which is a combination of 
the inversion after a data presmoothing and regularized collocation. In prin¬ 
ciple, each of the combined methods can be used alone as a regularization 
method. However, their performance very much depends on the interplay be¬ 
tween the degree of ill-posedness of Q, Q, coded in the rate of increase of 
the spherical symbol {a*,}, and the smoothness of the solution. 

For example, from Figure |T] one may conclude that for a severely ill-posed 
problem and moderately smooth solution (Figure [l£i) the regularized colloca¬ 
tion performs a bit better than the inversion after a data presmoothing, but 
for a moderately ill-posed problem (Figure [l]:) this is not the case. 

Since for ill-posed problems the smoothness of their solutions is usually 
unknown, it is not clear which of the single parameter regularizations is the 
better one. In this situation the combination of the single parameter regu¬ 
larizations, resulting in a two-parameter scheme, can be seen as a method of 
choice. The above figures illustrate this point. 

At the end of Section 4 we have used Theorems 2 and 3 to demonstrate 


theoretically that the two-parameter scheme (181 is able to provide a better 


error bound than ones for the schemes involved in the combination. However, 
this demonstration is limited to an a priori choice of the regularization param¬ 
eters, requiring knowledge of the solution smoothness, and therefore seldom 
usable in practice. 
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On the other hand, in Section 5 the above mentioned ability of the two- 
parameter scheme has been demonstrated for an a posteriori choice of the 
regularization parameters, where no knowledge of the solution smoothness has 
been used. At the same time, a posteriori parameter choice employed in Sec¬ 
tion 5 is of a heuristic nature. We hope that this paper will stimulate research 
towards theoretical justification of a posteriori rules for multi-parameter reg¬ 
ularization. 
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